10x-180831 Monocle analysis
library(Seurat)## Loading required package: ggplot2## Loading required package: cowplot##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggplot2':
##
## ggsave## Loading required package: Matrixlibrary(monocle)## Loading required package: Biobase## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:Matrix':
##
## colMeans, colSums, rowMeans, rowSums, which## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colMeans,
## colnames, colSums, do.call, duplicated, eval, evalq, Filter,
## Find, get, grep, grepl, intersect, is.unsorted, lapply,
## lengths, Map, mapply, match, mget, order, paste, pmax,
## pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
## rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which, which.max, which.min## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: VGAM## Loading required package: stats4## Loading required package: splines## Loading required package: DDRTree## Loading required package: irlbalibrary(pheatmap)
seurobj <- readRDS('../../10x-adipocyte-analysis/output/10x-180831')DE genes res.0.5 T1T2T3 and T4T5 combined - unregressed/unfiltered data
cds <- readRDS('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/monocle/180831/monocle_DEgenes-res.0.5-noreg-T1T2T3-T4T5-combined/10x-180831-noreg-monocle')plot_cell_trajectory(cds, color_by = 'timepoint')
DE genes res.1.5 T1T2T3 and T4T5 combined - unregressed/unfiltered data
cds <- readRDS('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/monocle/180831/monocle_DEgenes-res.1.5-noreg-T1T2T3-T4T5-combined/10x-180831-noreg-monocle')fig <- plot_grid(ncol=1,
plot_cell_trajectory(cds, color_by='timepoint'),
plot_cell_trajectory(cds, color_by='Pseudotime'),
plot_cell_trajectory(cds, color_by='State'),
plot_cell_trajectory(cds, color_by = "State") + scale_color_manual(values=c("#f67770", "#964B00", "orange"), name = "State"))
fig
#save_plot('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/monocle_plots.pdf', fig, base_width=5, base_height=16)
#save_plot('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/monocle_timepoint.png', plot_cell_trajectory(cds, color_by='timepoint'), base_width=7.5, base_height=6)
#save_plot('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/monocle_pseudotime.png', plot_cell_trajectory(cds, color_by='Pseudotime'), base_width=7.5, base_height=6)
#save_plot('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/monocle_state.png', plot_cell_trajectory(cds, color_by='State'), base_width=7.5, base_height=6)
#save_plot('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/monocle_prediction.png', plot_cell_trajectory(cds, color_by = "State") + scale_color_manual(values=c("#f67770", "#964B00", "orange"), name = "State"), base_width=7.5, base_height=6)plot_cell_trajectory(cds, color_by = "timepoint") + geom_point(color='white', size=5) + geom_point(aes(colour=timepoint), alpha=0.1)
plot_cell_trajectory(cds, color_by = "timepoint") + geom_point(color='white', size=5) + geom_point(aes(colour=timepoint), alpha=0.01)
#data <- readRDS('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/10x-180831-noreg')
#data <- AddMetaData(data, pData(cds)['State'])
#saveRDS(data, '/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/10x-180831-noreg')BEAM
BEAM takes as input a CellDataSet that’s been ordered with orderCells and the name of a branch point in the trajectory. It returns a table of significance scores for each gene. Genes that score significant are said to be branch-dependent in their expression.
#BEAM_res <- BEAM(cds, branch_point = 1, cores = 1)
load('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/monocle/180831/BEAM')
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]paste('Significant genes with q-val < 0.01:', length(BEAM_res$qval[BEAM_res$qval < 0.01]))## [1] "Significant genes with q-val < 0.01: 8360"paste('Significant genes with q-val < 0.001:', length(BEAM_res$qval[BEAM_res$qval < 0.001]))## [1] "Significant genes with q-val < 0.001: 7216"paste('Significant genes with q-val < 0.0001:', length(BEAM_res$qval[BEAM_res$qval < 0.0001]))## [1] "Significant genes with q-val < 0.0001: 6421"paste('Significant genes with q-val < 0.00001:', length(BEAM_res$qval[BEAM_res$qval < 0.00001]))## [1] "Significant genes with q-val < 0.00001: 5860"paste('Significant genes with q-val = 0:', length(BEAM_res$qval[BEAM_res$qval == 0]))## [1] "Significant genes with q-val = 0: 329"Histograms of p-values and q-values
hist(BEAM_res$pval)
hist(BEAM_res$qval)
You can visualize changes for all the genes that are significantly branch dependent using a special type of heatmap. This heatmap shows changes in both lineages at the same time. It also requires that you choose a branch point to inspect. Columns are points in pseudotime, rows are genes, and the beginning of pseudotime is in the middle of the heatmap. As you read from the middle of the heatmap to the right, you are following one lineage through pseudotime. As you read left, the other. The genes are clustered hierarchically, so you can visualize modules of genes that have similar lineage-dependent expression patterns.
# branched_8 <- plot_genes_branched_heatmap(cds[row.names(subset(BEAM_res,
# qval < 0.0001)),],
# branch_point = 1,
# num_clusters = 8,
# cores = 10,
# show_rownames = F,
# return_heatmap = T,
# branch_labels = c("Cell fate 1 (State 2)", "Cell fate 2 (State 3)"),
# branch_colors = c('#f67770', '#1bb840', '#649efc')
# )Heatmaps for all genes with q-val < 0.0001 (6421 genes), 6, 8 and 10 clusters.
load('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/branched2') Nr of genes per cluster
print('6 CLUSTERS')## [1] "6 CLUSTERS"for (i in 1:length(unique(branched_6$annotation_row$Cluster))){
cluster <- rownames(branched_6$annotation_row)[branched_6$annotation_row$Cluster == i]
print(paste('Nr of genes in cluster ', i, ': ', length(cluster), sep=''))
}## [1] "Nr of genes in cluster 1: 1510"
## [1] "Nr of genes in cluster 2: 900"
## [1] "Nr of genes in cluster 3: 1659"
## [1] "Nr of genes in cluster 4: 1375"
## [1] "Nr of genes in cluster 5: 544"
## [1] "Nr of genes in cluster 6: 433"print('8 CLUSTERS')## [1] "8 CLUSTERS"for (i in 1:length(unique(branched_8$annotation_row$Cluster))){
cluster <- rownames(branched_8$annotation_row)[branched_8$annotation_row$Cluster == i]
print(paste('Nr of genes in cluster ', i, ': ', length(cluster), sep=''))
}## [1] "Nr of genes in cluster 1: 1510"
## [1] "Nr of genes in cluster 2: 900"
## [1] "Nr of genes in cluster 3: 1117"
## [1] "Nr of genes in cluster 4: 1375"
## [1] "Nr of genes in cluster 5: 544"
## [1] "Nr of genes in cluster 6: 239"
## [1] "Nr of genes in cluster 7: 194"
## [1] "Nr of genes in cluster 8: 542"print('10 CLUSTERS')## [1] "10 CLUSTERS"for (i in 1:length(unique(branched_10$annotation_row$Cluster))){
cluster <- rownames(branched_10$annotation_row)[branched_10$annotation_row$Cluster == i]
print(paste('Nr of genes in cluster ', i, ': ', length(cluster), sep=''))
}## [1] "Nr of genes in cluster 1: 1510"
## [1] "Nr of genes in cluster 2: 900"
## [1] "Nr of genes in cluster 3: 1048"
## [1] "Nr of genes in cluster 4: 69"
## [1] "Nr of genes in cluster 5: 1375"
## [1] "Nr of genes in cluster 6: 161"
## [1] "Nr of genes in cluster 7: 239"
## [1] "Nr of genes in cluster 8: 194"
## [1] "Nr of genes in cluster 9: 542"
## [1] "Nr of genes in cluster 10: 383"# saveBEAMplot <- function(heatmap_object, filename){
# save_plot(paste("/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/BEAM_results/", filename, sep=''), heatmap_object, base_width=5, base_height=6)
# }
#
# saveBEAMplot(branched_6$ph_res, 'BEAM_qval0.0001_6clusters.png')
# saveBEAMplot(branched_6$ph_res, 'BEAM_qval0.0001_6clusters.pdf')
#
# saveBEAMplot(branched_8$ph_res, 'BEAM_qval0.0001_8clusters.png')
# saveBEAMplot(branched_8$ph_res, 'BEAM_qval0.0001_8clusters.pdf')
#
# saveBEAMplot(branched_10$ph_res, 'BEAM_qval0.0001_10clusters.png')
# saveBEAMplot(branched_10$ph_res, 'BEAM_qval0.0001_10clusters.pdf')6 clusters
branched_6$ph_res
8 clusters
branched_8$ph_res
10 clusters
branched_10$ph_res
Write BEAM results to files
writeBEAMtoFile <- function(branched_object, foldername){
for (i in 1:length(unique(branched_object$annotation_row$Cluster))){
BEAM_cluster <- BEAM_res[BEAM_res$gene_short_name %in% row.names(branched_object$annotation_row)[branched_object$annotation_row == i],]
write.table(BEAM_cluster, paste('/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/BEAM_results/', foldername, '/cluster', i, '.txt', sep=''), row.names=F, quote=F, sep='\t')
}
}
writeBEAMtoFile(branched_6, '6clusters')
writeBEAMtoFile(branched_8, '8clusters')
writeBEAMtoFile(branched_10, '10clusters')GSEA
Genes plotted over pseudotime
cds_subset <- cds[row.names(subset(fData(cds), gene_short_name %in% c("FABP5", "SCD", "G0S2", "ADIPOQ", "APOD", "MGP", "PLAC9", "DCN", "MT-CYB", 'CYB5A'))),]
#plot_genes_in_pseudotime(cds_subset, color_by="timepoint")
plot_genes_branched_pseudotime(cds_subset, branch_point = 1, color_by = "timepoint", ncol = 2)## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
save_plot("/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/genes_in_pseudotime.png", plot_genes_branched_pseudotime(cds_subset, branch_point = 1, color_by = "timepoint", ncol = 2), base_width=10, base_height=15)## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axisC/EBPa more important in white. C/EBPb and C/EBPd more important in brown (described in several reviews). Weird to see it the other way around in our data.
cds_subset <- cds[row.names(subset(fData(cds), gene_short_name %in% c('EBF2', 'PDGFRA', 'PDGFRB', 'PPARG', 'MALAT1', 'NEAT1', 'PRDM16', 'CEBPA', 'CEBPB', 'CEBPD', 'UCP1', 'LEP'))),]
#plot_genes_in_pseudotime(cds_subset, color_by="timepoint")
plot_genes_branched_pseudotime(cds_subset, branch_point = 1, color_by = "timepoint", ncol = 2)## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
save_plot("/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/genes_in_pseudotime.png", plot_genes_branched_pseudotime(cds_subset, branch_point = 1, color_by = "timepoint", ncol = 2), base_width=10, base_height=15)## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axiscds_subset <- cds[row.names(subset(fData(cds), gene_short_name %in% c('EBF2', 'PDGFRA', 'PDGFRB', 'PPARG', 'MALAT1', 'NEAT1', 'PRDM16', 'CEBPA', 'CEBPB', 'CEBPD', 'UCP1', 'LEP'))),]
plot_genes_branched_pseudotime(cds_subset, branch_point = 1, color_by = "State", ncol = 2)## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
save_plot("/projects/pytrik/sc_adipose/analyze_10x_fluidigm/data/plots_slides/genes_in_pseudotime2.png", plot_genes_branched_pseudotime(cds_subset, branch_point = 1, color_by = "State", ncol = 2), base_width=10, base_height=15)## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis